new
sel beam id=1 begin 0 0 0 end 0 0 10 nseg=20
sel beam prop emod=210e9 nu=0.3 xcarea=1.0 & 
xciy=0.0833 xciz=0.0833 xcj=0.167
sel node fix x y z xr yr zr range x=(-0.1,0.1) y=(-0.1,0.1) z=(-0.1,0.1)
sel beam apply zdist 1.0
sel set damp combined
pl sel beam moment my axes
solve

;---------Ԫӳ(ڵԪ)-----
def force_outp
jj=0
loop cid(1,20)
  sp=s_find(cid)
  if s_type(sp)=1 then
    jj=jj+1
    ;--1ڵ---
    np = s_node(sp,1)
    xx1 = nd_pos(np,2,1)
    yy1 = nd_pos(np,2,2)
    zz1 = nd_pos(np,2,3)
    ffx1= sb_force(sp,1,1) 
    ffy1= sb_force(sp,1,2) 
    ffz1= sb_force(sp,1,3) 
    mmx1= sb_mom(sp,1,1) 
    mmy1= sb_mom(sp,1,2) 
    mmz1= sb_mom(sp,1,3) 
    ;--2ڵ------
    np = s_node(sp,2)
    xx2 = nd_pos(np,2,1)
    yy2 = nd_pos(np,2,2)
    zz2 = nd_pos(np,2,3)
    ffx2= sb_force(sp,2,1) 
    ffy2= sb_force(sp,2,2) 
    ffz2= sb_force(sp,2,3) 
    mmx2= sb_mom(sp,2,1) 
    mmy2= sb_mom(sp,2,2) 
    mmz2= sb_mom(sp,2,3) 
    ; -----
    xx=(xx1+xx2)/2.0
    yy=(yy1+yy2)/2.0
    zz=(zz1+zz2)/2.0
    ffx=(ffx1+ffx2)/2.0
    ffy=(ffy1+ffy2)/2.0
    ffz=(ffz1+ffz2)/2.0
    mmx=(mmx1+mmx2)/2.0
    mmy=(mmy1+mmy2)/2.0
    mmz=(mmz1+mmz2)/2.0
    ;------------------------
    xtable(1,jj)=zz
    ytable(1,jj)=mmy
  endif
endloop
end

force_outp

set logfile beam_moment.txt
set log on 
set pagelength 10000
pri table 1
quit
